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f Abstract 

Oh 
.^^ We consider the interactions of finite dipoles in a doubly-periodic domain. A finite dipole is a pair of equal 

and opposite strength point vortices separated by a finite distance. The dynamics of multiple finite dipoles in an 
unbounded inviscid fluid was first proposed by Tchieu, Kanso & Newton in [l] as a model that captures the "far- 
field" hydrodynamic interactions in fish schools. In this paper, we formulate the equations of motion governing the 
dynamics of finite-dipoles in a doubly-periodic domain. We show that a single dipole in a doubly-periodic domain 
exhibits periodic and aperiodic behavior, in contrast to a single dipole in an unbounded domain. In the case of 
two dipoles in doubly-periodic domain, we identify a number of interesting trajectories including collision, collision 
O avoidance, and passive synchronization of the dipoles. We then examine two types of dipole lattices: rectangular and 

f~^ diamond. We verify that these lattices are in a state of relative equilibrium and show that the rectangular lattice is 

• unstable while the diamond lattice is linearly stable for a range of perturbations. We conclude by commenting on 

^ the insights these models provide in the context of fish schooling. 

>-. 

Ph 1 Introduction 

" ' The question of how interactions among individual fish result in highly-coordinated motion in fish schools has 

^^ been the focus of numerous studies, the majority of which consider behavior-based models of homogeneous 

J^ particles (fish) interacting locally based on rules of repulsion, alignment and attraction to other fish (see, for 

■<!;;j- example, |2] and (3). These models are capable of exhibiting realistic dynamics similar to those observed in 

O 

fT^ biological schools (see [1]), but do not elucidate the mechanisms by which individuals transmit and integrate 

• • information from the school to guide their motion as noted in [5 . In particular, little is known about the 

k>( role of the fluid medium in guiding the motion of the individual fish. Our main motivation in this paper is 

S to develop a framework for studying the hydrodynamic interactions inside a large school of fish, as opposed 

to near the school boundary, see Figure lllleft). We assume the school is homogeneous and we focus on fish 

interactions in a domain within the school. More specifically, we consider fish interactions in a rectangular 

domain with periodic boundary conditions. This argument holds when the characteristic length L of the 

domain is small relative to the school size but nuich larger than the fish size I and the separation distance 

R between two neighboring fish. 

As a leading-order model of the fish motion, we use N finite dipoles in a doubly-periodic domain, see 
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Figure 1: Schematic showing a homogeneous school of fish {left). To highlight fish interactions in a domain within the school 
as opposed to on its boundary, we use fish in a domain with periodic boundary conditions (middle). Individual fish are modeled 
using the finite dipole model (right) . 



Figure m middle). A finite dipole is a pair of equal and opposite strengtli point vortices separated by a finite 
distance, see Figure IIJ right), and thus its self-propelled speed is well-defined as opposed to that of a point 
dipole. For point dipole models, see [B], [7], [S] and [5]. Another reason for using the finite-dipole model is 
that a body propelling itself in a two-dimensional inviscid fluid produces, to leading order, a dipolar velocity 
field. For several hydro-dynaniically coupled swimming bodies, when their separation distance R is large 
relative to their size £, [Ij argued that a dynamical system consisting of a collection of interacting finite- 
sized self-propelled dipoles would be a reasonable model governing the far-field hydrodynamic interactions 
among the bodies. Numerical evidence in fB suggests that the finite dipole model is a good approximation 
of interacting bodies even when the separation distance i? ~ 3-^ is of order £. 

In this paper, we generalize the formulation of [1 to the case of a doubly-periodic domain. Given 
that each finite dipole consists of two constrained point vortices, we build upon known results on vortex 
interactions in doubly-periodic domains. The basic formulation of a simple vortex lattice, which is a special 
case of a vortex in a doubly-periodic domain, was first discussed by [10 and was extended by [IT] to general 
lattices. Vortices in periodic and doubly-periodic domains were further studied in [12], [13] and [M]. Clusters 
of point vortices in doubly-periodic domains were also considered in [15j . 

The organization of this paper is as follows. Section [2] addresses the formulation of N finite dipoles 
in a doubly-periodic domain as a constrained 2N point vortex system. This is done by directly modifying 
the standard point vortex equations of motion to respect the constraint that each pair of point vortices of 
equal and opposite strength remain a fixed distance £ apart. Section [3] discusses the periodic and aperiodic 



behavior of one dipole in a doubly-periodic domain. Examples of interactions of two finite dipole systems 
are presented in section [4J Section [5] focuses on two types of dipole lattices: rectangular and diamond. We 
show that both lattices are in a state of relative equilibrium with the former being unstable while the latter 
is linearly stable to a range of small perturbations. We conclude by commenting on the insights of these 
models offer in the context of schooling of fish. 

2 Problem Formulation 

Consider N pairs of point vortices or dipoles of equal and opposite strengths (±r„) placed a distance £„ 
apart, n = 1,. .. ,N. See Figure [l| middle) for a depiction of N dipoles in a doubly-periodic domain and 
Figure llTright) as well as for the details of a single finite dipole. The vortex of strength +r„ is referred to 
as the left vortex and its position is denoted by z„j whereas the — r„ vortex is called the right vortex and 
its position is denoted by Zn^,-. For convenience, complex notation {z — x + iy and i = \/— 1) is employed. 
The position z„ of the dipole center is related to z„j and z„_r via 

Zn = ^ • (1) 

Let a„ represents the orientation of the dipole with respect to the a;-axis. Then, the position of the left and 
right vortices is given by, 

Our goal is to formulate the equations of motion governing the interaction between N dipoles in a 
doubly-periodic domain while the vortex pair in each dipole is constrained to have a finite length £„ {£„ = 0) , 
hence the name finite dipole. This amounts to deriving equations of motion for all dipole centers z„ and 
the orientations «„. Following [Ij, we assume that the constraint £„ = induces an additional inter-dipole 
velocity to the 2Af- vortex problem. That is to say, we assume that the left and right vortices are advected 
according to the velocity 

Zn,i = Wn,s + Wn,o(2„j) + iA„e"'"" , (3a) 

Z„,r = Wn,s + Wn,oiz„,r) " iA„e"'"" , (3b) 

where ( ) and ( ) represent the complex conjugate and time derivative, respectively, and A„ is a real constant. 
The term Wn^s represents the sum of self-induced conjugate velocity of each dipole n and conjugate velocity 



induced by the dipole's own images. Note that due to the doubly-periodic nature of the domain, each dipole 
has infinitely many images. The term Wn^o represents the conjugate velocity induced by all other finite 
dipoles and their images. 

The term ±iA„e^'"" in ^ is an additional, attractive (A„ > 0) or repulsive (A„ < 0) inter-dipole 
velocity that allows us to apply the finite- length constraint on £„. Physically speaking, this term forces the 
two vortices z„.i and Zn^r to stay a fixed distance apart by allowing them to overcome the tendency to move 
toward or away from each other with speed |A„| along the line joining these two vortices. The introduction 
of non-zero A„ can be thought of as introducing a degree of freedom to enforce the constraint that £„ — 0, 
much like applying Lagrange multipliers in constrained mechanics. This allows each finite dipole to retain 
its 'particle-like' identity throughout its time evolution. It is noted in [T] that A„ translates into constraint 
forces acting on the Euler equation governing the fiuid system, thus breaking the Hamiltonian nature of the 
system. 

Equations (|3| can be rewritten as a system of equations governing the motion of the center 2;„ and the 
orientation a„ of each finite dipole. This is done exactly as in 1 for the unbounded plane except that Un^s 
and ujn.o have different expressions in the doubly-periodic domain as will be shown below. Namely, upon 
substituting ^ into ([s]), one gets 

^ , Wn,o{^n.l) + Wn,o{Zn.r) , ,. 

Zn=Wn,s^ ^ — , (4) 

Re [{Wn.o{Zn,r) ~ W„,o(z«.l)) e'""] _^ 

"n = ■ -. ■ , (5) 

and 

A„ = -Im [{Wn^o{Zn,\) - Wn,o{Zn.i-)) e'""] . (6) 

In order to close this system of equations, one needs to find expressions for the self-induced velocity Wn^s 
and the velocity Wn^o induced by other dipoles in a doubly-periodic domain. 

To obtain expressions for the terms Wn.s and Wn,o, first remember that the velocity field created 
by unconstrained vortices in a doubly-periodic domain is given in terms of the Weierstrass-C function (see, 
for example, [TU], [O], [H] and [TS]). This yields, upon straightforward manipulations, that the conjugate 
velocity field created by the unconstrained but paired 2N point vortices in a doubly-periodic domain is given 



by 



z = 7r^[C{z~ z„,;;iLJi,a;2) - C,(z ~ Zny,toi,L^2) 

ZTTl 






where C,{z; 101,102) is the Weierstrass C-function, 
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Op5 = 2pa;i + 2ga;2 , p,geZ-{0}. 



(7) 



(8) 



Here, p and g are signed integers. In ([7]), wi and UJ2 are the half-periods of the doubly-periodic domain, rji 
is the value of the Weierstrass C- function at the half-period cui , and A is the area of the rectangular domain 
and is given by the Legendre's relation, A = 2i{ujiuJ2 — ajiiuj2). Note that the third term in (It]) goes to zero 
for a square domain {uji is purely real and u)2 is purely imaginary, such that wi — \i-02\) )• 
By virtue of ([7]), one can readily verify that 'Wn,s takes the form 
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(9) 



whereas the term Wn^o representing the velocity induced by all other dipoles and their images is given by 



N 



3=^n 



o{z) = X! 9^ C,{z- ZjX,^\,^2) -C(2-2j,r;wi,a;2) 



(10) 



Awi ujij ■' A ■' 



Equations ^, ( 10 ) can now be substituted back into Q, (Is]) to get a closed system of 3N real equations (N 
complex -|- N real) governing the motion of N finite dipoles interacting in a doubly-periodic domain. The 
strength of the Lagrange multiplier A„ is obtained by substituting ^, (10) into (l6|. 

Two remarks are in order here. First, when the period of the system goes to infinity, wi, W2 — > 00, the 
self-induced velocity and the velocity induced by other dipoles are reduced to their counterparts in the case 



of an unbounded plane, namely, 



r„e- 
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(11) 



^-^ 1'K\ \Z — Zjl Z — Z 



J:l- 



Second, in the unconstrained interaction of N dipoles (or equivalently, 2N point vortices) in a doubly-periodic 
domain, the inter-dipole spacing £„ is not constant and the equations of motion for each vortex are given by 
substituting Q, (10 1 into (pi) with A„ = 0. The system of equations is then Hamiltonian and the total linear 
impulse is conserved. The Hamiltonian of a 2iV-vortex system subject to periodic boundary conditions can 
be written as 



2^ 2N 



^ = ^^^Y1 ^«^j { ^^ i'^^^" ~ ^j)i 

f TTUJl ril\ (Z„ - Zj)'^ 



n=l j=l 



(12) 
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\Aaji LOi 



^[{xn ~ Xjf + {yn - yjf]y 
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where n ^ j, a{z) is the Weierstrass sigma function and C(z) is the logarithmic derivative of (t{z). The total 
linear impulse of a 2 A^- vortex system "^j^i^jZj is conserved, whereas the angular impulse '^j^i^jZjZj is 
not conserved in periodic domains. This Hamiltonian structure is destroyed in the finite dipole system due 
to the constraint £„ = and the associated Lagrange multiplier A„ as noted above. 



3 Periodic and aperiodic behavior of single dipole 

We consider the seemingly simple case of a single dipole in a doubly-periodic domain. The dipole is only 
subject to its self-induced velocity and the velocity induced by its own images. That is to say, Wn,o and A„ 
are identically zero and equations ^, ^ take the form 
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(13) 



Note that the second term on the right-hand side of the first equation is identically zero in a square domain. 
This equation can be integrated in closed form. 
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t + z(0), a = a(0). 



(14) 



Here, z(0) and a(0) are the initial position and orientation of the dipole respectively. That is to say, a single 
dipole always moves in a straight line with its orientation angle a unchanged. However, the slope of the 
linear trajectory of the dipole center is not in the direction of the orientation angle a, see Figure[2] except for 
special initial conditions such as a(0) = fc7r/2, k arbitrary integer. The reason for this discrepancy between 
the slope of the dipole trajectory and its orientation is due to the periodic effect brought by the dipole images 
included in the ('-function. 

A single dipole affords two distinct types of dynamical behavior: aperiodic and periodic. For the 
former type, the single dipole traces out path which fills up the whole domain as depicted in Figure ^a). 
Here, the dipole never returns to the same location it visited thus the term aperiodic. Note that in this 
and the coming sections, we consider a doubly-periodic domain with half period wi = 5 and i>Ji — 5i. In all 
simulations, the strength F of the dipoles' vortices is set to unity and the finite dipole length i = l/27r. To 
emphasize the doubly-periodic nature of the domain, the dipole trajectory is plotted on a torus as done in 
Figure p[b). The torus is obtained by applying the linear transformation {x/ioi^y/ |w2|) — > [utV) £ [— tt, tt] 
and (wi, 1^2 1) — )■ (i?, r), where u, v are the angles and R, r are the radii defining the torus. More specifically, 
a point {X,Y,Z) on the torus is given by X = (i? + r cosw) cosu, Y — {R + rcosv)sa\u and Z — rsin?;. 
The ratio of r to i? is selected to be 1 to 2 for clarity of exposition. 

The aperiodic behavior in Figure [3] seems to be the generic behavior for arbitrary initial conditions. 
We then ask for what initial conditions (if any), the single dipole exhibits periodic solutions. That is to say, 
we look for solutions that satisfy the condition 



z(r) =^{Q) + 2puJi+2qu2, 



(15) 



where p and q are integers and T is the period of the motion. Using (14). the above equation amounts to 
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T = 2pLUi + 2ga;2, p,q £ Z. 



(16) 



In a square domain, wi = —iuj2 = i-^, ( 16 1 is satisfied only when the imaginary and real part on the left hand 



side of the equation is in rational ratio of q to p. The corresponding value of a can be evaluated numerically 
for different q to p ratio. See Figure [4] for a depiction of a periodic trajectory for q/p = —2. 

We close this section by noting that in the limiting case of an infinitely large domain, wi, uj2 — > c». 



C(z) reduces to l/z and A— > in (15). Left hand side of (151 thus becomes 
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(17) 
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Figure 2: Comparison between trajectories of a single dipole in unbounded plane and in doubly-periodic domain for wi = 0.2, 
(jJ2 = 0.2i and o(0) = 7r/3. Black color denotes the path taken by the dipole in a doubly-periodic domain while grey color 
denotes the path taken by the dipole in unbounded plane. In both cases, the orientation a remains constant for all time but in 
the doubly-periodic domain, the slope of the trajectory traced by the dipole center is not equal to a. 
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Figure 3: Aperiodic trajectory that densely fills the whole domain: (a) trajectory depicted in doubly-periodic domain (b) same 
trajectory depicted on a torus. Parameter values are: i = 1/2it, uj = 5, z{0) = 0, ce{0) = tt/3. 



The ratio of the imaginary and real part of the above expression is -tan a and the self-induced velocity 



reduces to the case of an unbound plane, represented by ( 11 ). 



4 Collision, no-collision and synchronization of two dipoles 

In this section, we consider the interaction of two finite dipoles (zi, ai) and (z2, 02) of equal length £1 = ^2 == ^ 
and equal strength Fi = r2 = F, we describe three distinct dynamical behavior: collision, collision-avoidance 
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(a) Rational ratios of q/p give o(0) that produce periodic trajectories 
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Figure 4: Periodic trajectory of a single dipole in doubly-periodic domain: (a) according to equation \16\ , rational values of 
q/p give o(0) that produce periodic trajectories. (6) periodic trajectory for p/q = —2 and parameter values I = l/27r and ui = 5. 
The value of a(0) is obtained from plot (a). (6) same trajectory depicted on a torus. 



and motion synchronization of the two dipoles. These nontrivial interactions arise solely from hydrodynamic 
coupling. 

For concreteness, we write the equations of motion for the two dipole system by substituting (|9| and 
p^ into Q and ([5|. To this end, one has 



zi 



27ri 



{C(-i^c'"0 + 2 [C(^l,l - ^2,l) - C(2l,l - 22,r) + C(zi,r - ^2.l) - C(zi,r " Z2,,) 



i£ e 
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AuJl UJl 



i£(e-i"i+e^'"^)^}, 



(18) 



and 



3^1 = Rej— e'"i C(^l,r - Z2,\) - C(^l,r " 2:2,r) ^ C(^l,l - ^2,l) + ((^^l,! - ^2,r) | 



(19) 



Similar equations hold for Z2 and 012. These equations form a system of six coupled nonlinear ordinary 



differential equations which we solve numerically using a standard Runge-Kutta solver with variable time 
step. For a fixed set of parameter values ii ^ £2 ^ l/27r, wi — \uj2\ = 5, we vary the initial conditions 
zi(0), ai(0), 22(0) and 0:2(0). Depending on the choice of initial conditions, one obtains different dynamical 
behavior. Note that the periodicity of the domain enriches the dynamics of the two dipoles and enables them 
to interact many times: when one dipole leaves the domain, it re-enters on the opposite side and continues 
to interact with the other dipole. In this sense, the two dipoles cannot diverge but do exhibit a range of 
dynamical behavior as summarized below. 



y 




Figure 5: Collision of two dipoles in a doubly-periodic domain. The blue and red line represents the trajectories formed by 
the two dipoles. Parameter values are: i = 1/2tt, 21 (0) = —2.5, 22(0) = 1.5 — 2i, ai(0) = tt/S, 02(0) = 7r/12. 



Collision. Two finite dipoles in doubly-periodic domain may collide in finite time. When the dipoles 
collide, they form a fixed quadrupole. A typical example of dipole collision in a periodic domain is shown 
in Figure pa. A similar behavior is reported in ^ for two finite dipoles in an unbounded plane. However, 
the interactions of two dipoles in an unbounded plane is simpler in the sense that one can identify the set 
of initial conditions that give rise to collision. In the doubly-periodic domain, the set of initial conditions 
that lead to collision seems to be dense in the space of all initial conditions (based on a range of numerical 
simulations not shown here for brevity). 

Collision avoidance. The hydrodynamic coupling between two dipoles could induce collision avoidance as 
shown in Figure [6^. Figure [6|d plots the change in the dipoles' orientation as a function of time. As the two 
dipoles approach each other, their orientations change drastically and collision is avoided. The trajectories of 
the dipoles are represented on a torus in Figure [6b. It is important to emphasize that this collision avoidance 
behavior is a result of the hydrodynamic coupling only with no external control. That is to say, the fluid 
medium plays the role of a collision avoidance mechanism for certain approach conditions. Again, numerical 
evidence (results not shown here) suggests that the set of initial conditions leading to collision avoidance is 

10 



2.5 



y Ox 



-2.5 




-5 -2.5 2.5 5 

X 

(a) in doubly-periodic domain 




(b) change in angles of orientation with time 




(c) on torus 

Figure 6: Collision-avoidance of two dipoles. (a) trajectories of the two dipoles. The blue and red dashed line are the paths of 
the dipoles if they act independently without interacting with each other, (b) Orientation angles versus time, (c) trajectories 
depicted on a torus. Parameter values are: t = l/27r, ^i(O) = —4.5, .22(0) = —3.5 — 2.5i, «i(0) = 0, 012(0) = 7r/6. 

dense in the space of all initial conditions. 

Synchronization. The most remarkable interaction mode of the two dipoles is the synchronization mode. 
We use the term synchronization to denote periodic trajectories where the two dipoles oscillate and interact 
with each other and return to the same position periodically. These periodic trajectories are found numeri- 
cally using an iterative method. In particular, we use a shooting method that adjusts the initial position of 
one of the dipoles iteratively to hone in on the periodic orbit. 

We distinguish two types of synchronization trajectories: unbounded and bounded. For the unbounded 
mode, the dipoles move side-by-side along undulating paths that exit the doubly-periodic domain to re-enter 
on the other side as depicted in Figure [7J In the bounded mode, the two dipoles dance around each other 
tracing out flower like orbits as depicted in Figure [8| These dancing trajectories are bounded in the sense 
that the dipoles move within a confined region of the doubly-periodic domain. 
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Figure 7: (a) Periodic trajectories of two dipoles. 2i(0) = 0, ^2(0) = 0.1501 + 0.4901i, oi(0) = 0, 02(0) = 7r/12. (6) Same 
synchronization trajectories depicted on a torus. 

5 Stability of Rectangular and Diamond Lattices 

We address the dynamics and stability of two families of dipole lattices: rectangular and diamond (see 
Figure ^ . By lattice, we mean an arrangement of dipoles in an ordered pattern extending to infinity in 
the unbounded plane. A rectangular lattice consists of dipoles aligned along ak^m — '^/'^ with their centers 
placed at Zk^m = ka + imh, where k, m = 0, ±1, ±2, ±3, . . ., a denotes the distance between two neighboring 
dipoles of the same row, and b denotes the distance between two rows. A diamond lattice consists of 
dipoles aligned along ak^m — 7r/2 but with centers placed such that Zk^m — ka + imb, for k,m even, and 
Zk,m ^ {k + -)a + 1(771 + -)b, for k,m odd. An alternative, and perhaps more elegant way, of describing 
these lattices is by considering them as special cases of dipoles in doubly-periodic domains. We adopt the 
latter view in this section. In particular, we define the 'smallest' doubly-periodic domain (or 'smallest cell') 
needed to describe these rectangular and diamond lattices. We then use the formulation in sections [2] and [3] 
to prove that these configurations correspond to relative equilibria of the finite-dipole dynamical system and 
we analyze their linear stability. 

The 'smallest' doubly-periodic domain needed to generate the rectangular lattice has half-periods wi = 
a/2 and a;2 = ifo/2 and contains a single dipole with orientation q;(0) = 7r/2, see Figure l9](a). The dipole's 



center is placed at the center of the domain for convenience. One can readily verify, using equation (13 1, 



that, for all time, the dipole's orientation a remains unchanged while its center moves with constant velocity, 
thus the rectangular lattice is in a state of relative equilibrium. 

To generate the diamond lattice, the 'smallest' doubly-periodic domain has half-periods uji — a/2 and 
UJ2 = 16/2 and contains two finite dipoles (zi, ai) and (2:2, 02) of equal strength F and equal length ii — i2 = i, 
with orientations ai(0) = 012(0) = 7r/2 and positions ^2(0) — Zi{0) + (a/2 + 16/2), see Figure loFb). To prove 
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Figure 8: Dancing of two dipoles. In all cases, ^i(O) = 0, 01 (0) = and 012(0) = 57r/6. Value of Z2{0) in each case: (a) 
0.01952 + 0.2i, (b) 0.11 + 0.13007i, (c) 0.08 + 0.13989i, (d) 0.11 + 0.116011, (e) 0.04999 + 0.20091, (/) 0.01002 + 0.21. Arrows 
are only drawn in (a) and (b). The dipoles are tracing these trajectories in anti- clockwise direction. 
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that this configuration is a relative equihbrium of the equations of motion, one needs to show that, for all 
time, zi = 2,2 is constant and 6i\ — a-i — 0. Actually, it suffices to show that Z2 — zi = and ai = 0:2 = 



for all time. Given these conditions, it immediately follows from ( 18 1 that the lattice's velocity z\ — 2,2 is 
constant for all time. The fact that the relative velocity Z2 — zi = is zero is obtained using ( [T8| (and its 
analog for Z2) to get 



Z2 - Zi 
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C(-Z2,l - 2:1,1) - C(-Z2,l - Zl,r) + C('Z2,r - ^l.l) " C(22,r - ^l.r) 

- C(2l,l - -Z2,l) + ((^1,1 - 22,r) - C(2l,r - 22,l) + CC^l.r ^ 22,r) 



(20) 



Now, recall that the C-function is an odd function (that is to say, C,{z) = — C(— z)) and note that, for the 



diamond lattice, one has Z2 — zi ~ Z2_; — zu — Z2,r ~ -^i.r = ^i +1^2- Then, it follows immediately from (20) 



that Z2 — zi ^ 0. Similarly, to show that the dipoles' orientation is constant for all time, that is to say, that 



di = 0:2 = 0, rewrite ( 19 1 in the form 



di = Re 
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Z1-Z2- ^(e-^ + e-^)) + azi -Z2 + ^(e'"^ + e'"^)) - 2C(zi - Z2; 



(21) 



Now, recall the periodicity property of the ^-function and the fact that, for the diamond lattice zi — Z2 
—uJi — UJ2, one can readily verify that the following identities hold 



C(zi - Z2) = c (-^1 - ^2) = -VI - m, C{zi -Z2 + ife'") = C(-wi - ^2 + ife'"), 



C(zi - Z2 - ife'") = -C{-^i - CJ2 + ifc'") - 2-qi - 2%. 



(22) 



Substitute (22) into (21) and use the fact that ai(0) = a2(0) and zi — Z2 is constant to get that di(i) = 
di(0) = 0. The same result holds for d2. Thus, the diamond lattice is a relative equilibrium of the finite- 
dipole dynamical system. 

Let C/iattico be the constant translational velocity of the lattice. It is instructive to compare the lattice 
velocity of the diamond and rectangular configurations with the velocity of a single dipole in an unbounded 
domain. For F = 1 and I = l/27r, one can readily see using (11) that the velocity of a single dipole in an 



unbounded domain is equal to 1. We use the same parameter values and compute C/iatticc from ( 13 ) and ( 18 ) 



for the rectangular and diamond lattices, respectively. The results are shown in FigurefTolas a function of the 
lattice density, defined as the number of dipoles per unit square length. Obviously, neither the rectangular 
nor the diamond formation present any advantages over the single dipole in term of increased translational 
velocity. To the contrary, the hydrodynamic interactions cause the dipoles in these lattices to move slower 
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Figure 9: Schematic of (a) rectangular lattice and (6) diamond lattice. The smallest doubly-periodic domain that generates the 
lattice is depicted in light grey. A larger doubly-periodic domain or "cell" is also depicted in solid red. 

than the single dipole with the rectangular lattice being slowest for all shown densities. A further increase 
in the lattice density causes C/iatticc to reverse sign and the dipoles to move in the direction opposite to their 
self-induced velocity. By continuity arguments, one deduces that there exist critical density values for which 
the rectangular and diamond lattices are stationary. 

We now examine the linear stability of these relative equilibria. Typically, the stability of infinite 
lattices is analyzed by introducing infinitesimal perturbations on each dipole's position and orientation, and 
looking for plane wave solutions of the linearized system; see, for example, 1 16] for a review of the stability 
analysis of a row of point vortices and of a von Karman street. See also fT7\ for stability of 2D vortex 
lattices and the more recent work [18 on the stability of driven and motile particle lattices in confined 
geometry. This approach involves infinite sums whose convergence needs to be established. Here, we avoided 
this complication by using a doubly-periodic domain and the Weierstrass C-function. Indeed, in [19], Aref 
showed that the stability of an infinite row of point vortices can be formulated and studied as the stability 
of point vortices in a periodic domain. He noted that the perturbation wave solution is equivalent to 
the eigenvalue problem associated with point vortices in a periodic domain. Further, since a wave of any 
wavelength must repeat after a finite number of vortices, various wavelengths can be captured by considering 
vortices in a periodic domain of various "cell" sizes. A cell is a doubly-periodic domain that is not necessarily 
the smallest, as depicted in Figure [9j We follow Aref's approach in the sense that we consider dipoles in a 
doubly-periodic domain and we apply perturbations in cells of various sizes to analyze how the stability of 
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Figure 10; Lattice translational velocity (/lattice versus its density for a square doubly-periodic domain \uji\ = |aj2| = io, and 
parameter values F = 1, ^ = 1/2it. Square symbol corresponds to the rectangular lattice while diamond symbol corresponds to 
the diamond lattice. A single dipole in an unbounded domain has unit velocity, shown in straight solid line. 



the lattice depends on the periodicity of the perturbation. 

For concreteness, we consider a doubly-periodic cell containing N dipoles. Let Sz^ and Sa^, n = 
1,2, . . . , N , denote the infinitesimal perturbations on the position and orientation of each dipole in this 



[^lattice + ^z„ and a„ = f (5a,i. Due to the doubly-periodic nature of the problem, 



cell so that z„ 

these perturbations will be repeated periodically. We linearize equations Q and ([5]) about the unperturbed 
lattice configuration and make use of the formula dC,{z)/ dz ~ —p{z), with p(z) being the Weierstrass Elliptic 
function defined as 



p{z;uji,uJ2) 



1 



E 



1 



yZ liLpqj 



1 

"P9 



p,qeZ-{0}. 



The linearized perturbed equations can be written in matrix form as follows 



di 



( X \ 

OXi 




bXj 


Syi 


= M^J 


^Vj 


[Sa, ) 




\ Sa, I 



hJ 



,iV 



(23) 



(24) 



The eigenvalues of the Mij matrix are computed numerically for all i,j. The lattice is said to be linearly 
stable to a given perturbation if all eigenvalues of Alij have non-positive real parts, that is to say, if all 
Re(A)< 0. The analysis is performed systematically by considering doubly-periodic cells of various sizes, 
starting with the smallest domain size a, b. When the cell size is equal to the smallest domain, the same 
perturbation is applied to all dipoles. The values of the largest Re(A) are tabulated in Table 1 for the 
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Table 1: Rectangular lattice. 

largest Re(A) 





a= 1 


6 = 


1 


Cell Size 


N 


6=1 


6= 1.5 


6 = 2 


a = 1.5 


a = 2 


a,b 


1 

















2a, 2b 


4 


1.09 


0.84 


0.80 


0.84 


0.80 


3a, 36 


9 


0.91 


0.73 


0.71 


0.73 


0.70 


4a, 46 


16 


1.09 


0.84 


0.80 


0.84 


0.80 


5a, 56 


25 


1.02 


0.80 


0.77 


0.80 


0.76 


6a, 66 


36 


1.09 


0.84 


0.80 


0.84 


0.80 



Table 2: Diamond lattice {b = 1) 





largest Re(A) 


Cell Size 


N 


a = 1 


a = 1.1 


a = 1.2 


a = 


= 1.3 a 


= 1.4 


a = 1.5 


a = 2 


a, 6 


2 

























2a, 26 


8 






















0.67 


3a, 36 


18 



















0.15 


0.58 


4a, 46 


32 



















0.15 


0.67 


5a, 56 


50 
















0.07 


0.13 


0.64 


6a, 66 


72 
















0.09 


0.15 


0.67 



rectangular lattice and Tables 2 and 3 for the diamond lattice. 

This analysis shows that the rectangular lattice is always unstable while the diamond lattice can 
be either unstable or linearly stable, depending on the lattice parameters a and 6 and on the size of the 
cell where the perturbation is applied. For example, when the same perturbation is applied to all dipoles, 
that is to say, when the size of the cell where the perturbation is applied is the same as the size of the 
smallest doubly-periodic domain, both the rectangular and diamond lattices are linearly stable. Also, when 
a = 6 = 1, the diamond lattice is always linearly stable but not the rectangular lattice. For a = 2, 6 = 1, the 
diamond lattice is always unstable (except as we just noted when the perturbation domain is the smallest 
doubly-periodic domain) . 

These results have been confirmed by numerically integrating the nonlinear equations in Q and ([5]) 
for the perturbed lattices. The perturbations are chosen randomly such that their magnitude is of the 
order a/1000. For the cases predicted to be unstable by the eigenvalue analysis, the lattices break down in 
finite time. Figure [TT] provides snapshots of the collapse of a rectangular lattice subject to initial random 
perturbations applied in the shown domain {N — 16). Meanwhile, for the linearly stable cases, the lattice 
keeps its integrity as shown in Figure [T2| for a diamond lattice with parameters a = 6 = 1. The diamond 
formation persisted to the end of the integration time (T = 100 time units). 

To quantify the deviation from the unperturbed lattice structure, we compare the dipoles positions 
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Table 3: Diamond lattice {a = 1) 







lai 


■gest Re(A) 






Cell Size 


N 


6 ==1.1 


6=1.2 


6 


= 1.3 6 


= 1.4 


6= 1.5 


a,b 


2 



















2a, 26 


8 













0.16 


0.42 


3a, 36 


18 













0.22 


0.37 


4a, 46 


32 













0.30 


0.42 


5a, 56 


50 










0.18 


0.33 


0.42 


6a, 66 


72 










0.20 


0.34 


0.45 



at each time t with that of the unperturbed lattices using 



e(i) 



Ep,, \zp{t)'Z,{t)f~\z]; 



lattice ^lattice I 

Zq I 



{N- l){N-2)/2 



PT^q- 



(25) 



This expression can be thought of as the mean square deviation of the perturbed lattice compared to the 
unperturbed lattice and its value is shown in Figure [13] for the perturbed rectangular and diamond lattices 



of Figures 11 and[T2j respectively. Clearly, the deviation of the perturbed rectangular lattice begin to grow 
rapidly around t — 3 time units, while the deviation of the perturbed diamond lattice remains small for all 
integration time. 

We conclude this section by commenting on the insights these results provide in the context of fish 
schooling. For fish schools, it has been argued that the diamond formation is favorable from an energy 
efficiency standpoint, |20| . In this seminal work, Weihs based his analysis on a stationary infinite diamond 
lattice and computed the locomotory benefits a given fish gets from the vortical wakes of neighboring fish. 
The wakes were modeled as idealized vortex streets and the fish were assumed to be point particles. That is 
to say, Weihs' model accounted for the near-field effect of fish wakes. The finite-dipole model considers the 
far-field hydrodynamic coupling (neglecting near-field vorticity) of self-propelled swimmers, and, as such, 
can be viewed as complementary to Weihs' model with the important difference that it allows for dynamic 
interactions among the fish (dipoles) whereas the latter assumes stationary fish. Based on the finite-dipole 
model, we make the following observations: 

(i) Neither the rectangular nor the diamond dipole lattices provide locomotory advantages to the individual 
dipoles in the sense that the lattice translational velocity is smaller than the velocity of a self-propelled 



dipole in an unbounded plane (see Figure 10). Perhaps not surprisingly, this result emphasizes that any 



locomotory advantages to schooling in terms of efficiency of motion would arise from near-field vortical 
wakes, as in Weihs' model, and not from far-field effects. Extraction of energy from near-field vorticity 
has been confirmed experimentally in live and dead trout, see pTll^ . 
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Figure 11: Time evolution of rectangular lattice for parameter values a = b= l,r = l and £ = 1/2tt. (a) initial configuration 
of rectangular lattice subject to random perturbations in the shown cell (N = IG). (6) and (c) trajectories of the dipoles at two 
different times, (d) Collision of dipoles and break down of the dipole lattice. 

(ii) However, our model shows that the diamond formation is beneficial from a stability standpoint. It is 
not clear how much stability is a desirable feature in large fish schools on the move. Stability, which 
measures how much a system opposes change, limits maneuverability. Here, we are referring to the 
stability and maneuverability of the school as opposed to that of the individual fish, the latter has been 
the topic of several studies, see, for example, [23] and references therein. We conjecture that passive 
stability of the school when subject to small perturbations might be desirable to a migrating school 
of fish. Active stabilization to stay in a school is energetically costly and therefore it may be more 
beneficial to travel in a school formation that is passively stable and requires no or little additional 
effort to maintain. These statements are yet to be validated by experimental observations. If true, they 
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Figure 12: Time evolution of diamond lattice for parameter values a = b=l, r = l and i = 1/2tt. (a) initial configuration of 
diamond lattice subject to random perturbations in the shown cell (N = 32j. (b) trajectories of the dipoles after an integration 
time T = 30 time units. 



£ 1 




Figure 13: t versus time t where e is the deviation of the perturbed lattices in Figures\ ll\and\l S\from their respective unperturbed 
structure. Clearly, the rectangular lattice looses its lattice structure whereas the diamond lattice maintains its lattice integrity. 

imply that the diamond formation is beneficial for both energy extraction from near-field wakes, |20) . 
as well as for passive stabilization of the school formation when subject to small perturbations. 



6 Conclusions 

We derived equations of motion for a system of finite dipoles in a doubly-periodic domain. We started from 
the standard point vortex equations in doubly-periodic domains and followed an approach similar to that 
in [T] for finite dipoles in unbounded plane. We used the resulting equations of motion to examine the motion 
of one and two dipoles in doubly-periodic domains. We showed that a single dipole in a doubly-periodic 
domain can exhibit periodic and aperiodic motion, whereas two dipoles exhibit a range of interesting behavior 
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including collision, collision-avoidance, and motion synchronization. In the latter category, the two dipoles 
travel in synchrony along unbounded and bounded periodic trajectories due to hydrodynamic coupling only. 
In the context of fish schooling, our main motivation for considering this class of models, these trajectories 
imply that hydrodynamic interactions may be responsible, at least in part, for the remarkable synchrony of 
motion observed in schools of fish. Further, the bounded periodic trajectories reported here are reminiscent 
to the stable epicyclic orbits that were observed in the context of vortex dipoles in the dilute-gas regime of 
a Bose-Einstein condensate, [24'. Indeed, it is known that equations governing quantized vortices in helium 
II and Bose-Einstein condensates are the same as that in ideal, incompressible fluids (see |1S|)- ^ formal 
connection between quantized vortices and the finitc-dipole model is beyond the scope of the present paper. 

We then identified two families of relative equilibria consisting of rectangular and diamond lattices, 
respectively. We examined the linear stability of these dipole lattices and found that the rectangular dipole 
lattice is always unstable whereas the diamond lattice is linearly stable for a range of parameter values 
and perturbation domains. In the context of fish schools, we argued that active stabilization to stay in the 
school formation is energetically costly and therefore it may be more beneficial for a migrating school of 
fish to travel in a diamond formation that is passively stable to small perturbations and requires no or little 
additional effort to maintain. 

Finally, we note that, motivated by recent advances in microfiuidics, the motion of self-propelled 
particles (bacteria) in two-dimensional fiuid channels (Hele-Shaw cells) has been the topic of several recent 
studies; see, for example, |18j and references therein. It is well-known that the equations of motion governing 
Hele-Shaw flows, though viscosity-driven, are identical to those of the inviscid potential flow. Therefore, we 
expect our model to be applicable in the microfluidic context as well. This direction will be pursued in future 
work. 
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